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Inverse transport problems in quantitative PAT for 

molecular imaging 

Kui Ren* Rongting Zhang^ Yimin Zhong^ 


Abstract 

Fluorescence photoacoustic tomography (fPAT) is a molecular imaging modality 
that combines photoacoustic tomography (PAT) with fluorescence imaging to obtain 
high-resolution imaging of fluorescence distributions inside heterogeneous media. The 
objective of this work is to study inverse problems in the quantitative step of fPAT 
where we intend to reconstruct physical coefficients in a coupled system of radiative 
transport equations using internal data recovered from ultrasound measurements. We 
derive uniqueness and stability results on the inverse problems and develop some effi¬ 
cient algorithms for image reconstructions. Numerical simulations based on synthetic 
data are presented to validate the theoretical analysis. The results we present here 
complement these in [57] on the same problem but in the diffusive regime. 

Key words. Photoacoustic tomography (PAT), molecular imaging, fluorescence optical tomog¬ 
raphy, fluorescence PAT (fPAT), radiative transport equation, hybrid inverse problems, numerical 
reconstruction 


1 Introduction 

Photoacoustic tomography (PAT) [10, 15, 20, 39, 41, 43, 47, 59, 66, 67, 68] is a recent hybrid 
imaging modality that attempts to reconstruct high-resolution images of optical properties of 
heterogeneous media. In a PAT experiment, we send a short pulse of near-infra-red (NIR) 
photons into an optically heterogeneous medium. The photons travel inside the medium 
following a radiative transport process. The medium absorbs a portion of the photons 
during their propagation process. The absorbed photons lead to the heating of the medium 
which then results in a local temperature rise. The medium expanses due to the temperature 
rise and then contracts when the rest of the photons leave the medium and the temperature 
drops accordingly. The expansion and contraction of the medium induces pressure changes 
which then propagate in the form of ultrasound waves. We then measure the ultrasound 
signals on the surface of the medium and from these measurements we intend to infer as 
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much knowledge as possible on the optical properties, for instance the optical absorption 
and scattering coefficients, of the medium. 

In recent years, there are great interests in developing PAT for biomedical molecular 
imaging [16, 51, 52, 65, 69, 68, 72], The main objective here is to visualize particular cellular 
functions and molecular processes inside biological tissues by using target-specific exogenous 
contrasts. To be specific, we consider in this work quantitative PAT for fiuorescence opti¬ 
cal imaging where one aims to image distribution of fluorescent biochemical markers inside 
heterogeneous media. In a typical imaging process, we first inject fluorescent markers into 
the medium to be probed. The markers will travel inside the medium and accumulate on 
their targets, for instance cancerous tissues inside the normal tissue. We then send a short 
pulse of NIR photons at wavelength Aa, to the medium to excite the fluorescent markers who 
then emit NIR photons at a different wavelength A^. The absorption of both the excita¬ 
tion and the emission photons by the medium will then generate ultrasound waves inside 
the medium following the photoacoustic effect just as in a regular PAT process, assuming 
that fluorescence takes place instantaneously as excitation light pulse is absorbed [57]. We 
then measure the ultrasound signals on the surface of the medium and attempt to recover 
information associated with the biochemical markers. 

The density distributions for the external light source and the fluorescent light in the 
tissues are both described by the radiative transport equation. Let Q {d > 2) be the 

domain of interests and be the unit sphere in We denote by X = 12 x the 

phase space and r± = {(x, v) G dfl x ± n(x) • v > 0} its boundary sets. We denote 

by Ua;(x, v) and Umi'x., v) the density of photons at the excitation and emission wavelengths 
respectively, at location x, traveling in direction v G Then Ua;(x, v) and Um(x, v) solve 

the following coupled system of radiative transport equations 

V ■Vu^ + {aa,^ + as^^)u^ = as,xKe{u^), in A 

V ■VUm + icra,m +Crs,m)Um = (T*,™ AT© (u^) TycTa,,,/ (x) A/(x), ffiX (1) 

Ux{x,v) = gx{x,v), u^(x,v) = 0, on r_ 

where the subscripts x and m denote the quantities at the excitation and the emission 
wavelengths, respectively. The coefficients aa^x and ag^x (resp. aa^m and ag^m) are respectively 
the absorption and scattering coefficients at wavelength Aa, (resp. A^). The scattering 
operator Kq and the averaging operator K/ are defined respectively as 

Ke(ux)(x,v) = / G(v,v')ux(x,v')dv', and, Kf(ux)(x,v) = / Ux(x,v')dv', (2) 

Jsd-1 

with the scattering kernel 0(v,v') describing the probability that a photon traveling in 
direction v' gets scattered into direction v. 

The total absorption coefficient aa,x consists of a contribution aa,xi from the intrin¬ 
sic tissue chromophores and a contribution (la^xf from the fluorophores of the biochemical 
markers: aa,x = cra,xi + cra,xf- The absorption coefficient due to fluorophores, (la^xf is propor¬ 
tional to the concentration p(x) and the extinction coefficient £(x) of the fluorophores, i.e. 
(^a,xf = ^(x)p(x). The coefficient ry(x) is the quantum efficiency of the fluorophores. The 
coefficients ry and (la^xf are the main quantities associated with the biochemical markers. 

The energy absorbed by the medium and the markers consists of two parts. The first 
part is from the excitation photons. This part can be written as <Ja,xKi{ux). The second part 
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of absorbed energy comes from emission photons. This part can be written as 

Therefore, the pressure held generated by the photoacoustic effect can therefore be written 

as: 


i 7 (x) = E:(x) (cr„,3,(x) - r7(x)cr„,3,/(x))/f/(n3,)(x) + cr„,^(x)/f/(n^)(x) 



( 3 ) 


where 5 is the (nondimensional) Griineisen coefficient that measures the photoacoustic 
efficiency of the underlying medium, and is the short notation for aa^^i + (1 “ v)^a,xf- 
We want to emphasize that when calculating the initial pressure held generated, we have 
subtract a portion of the energy, rjaa^xfKiiux), from the total energy absorbed by the medium 
and the markers. This is because that portion of energy is used to generate huorescence, 
not the heating in the photoacoustic process. 

The initial pressure held generated from the photoacoustic effect, H, evolves in space 
and time following the acoustic wave equation [14, 25, 62]: 



( 4 ) 


where c(x) is the speed of the ultrasound in the medium. The data that we measure are 
the solutions to the wave equation (4) on the surface of the medium, P|(o,tmax)x 9 Q, ^max being 
large enough, for various excitation light sources. 

Following [57], we call the process of reconstructing information on r] and (Ta,xf from 
datum P\(o,tmax)xdn huorescence PAT (fPAT). This is a molecular imaging modality that 
combines PAT with huorescence optical imaging. We refer interested readers to [57] for 
more discussions on the mathematical modeling of fPAT, including detailed derivation and 
justihcation the models (1) (in diffusive regime) and (4), and to [16, 51, 52, 65, 69] for some 
experimental and computational results on fPAT. Recent progress on huorescence optical 
imaging itself can be found in [5, 8, 27, 42, 48, 60] and references therein. 

Image reconstruction in fPAT is a two-step process as in regular PAT. In the hrst step, 
we reconstruct H from measured acoustic data. We assume here that this step has been 
hnished with methods such as those in [4, 7, 17, 18, 24, 30, 31, 34, 37, 40, 50, 62] and we are 
given the internal datum (3). Moreover, we assume that: (A-i) the Gruneisen coefficient 
5 as well as the absorption and scattering coefficients of the medium at the excitation 
wavelength, aa^^i and have been known from other imaging technologies (for instance 
a multi-spectral quantitative PAT step [13, 44]) before the huorescent biochemical markers 
are injected into the medium; and (A-ii) the absorption and scattering coefficients at the 
emission wavelength, aa,m and ag^mi are also reconstructed by other imaging methods (for 
instance a regular quantitative PA.T technique [6, 11, 12, 13, 14, 19, 21, 26, 44, 49, 56, 58, 73] 
after the Gruneisen coefficient is known). Therefore, our main objective is only to reconstruct 
the quantum efficiency r} and the fluorescence absorption coefficient aa,xf{^) in the system (1) 
from the datum H in (3). This is the quantitative fPAT (QfPAT) problem. 

Let us now remark on a couple of issues regarding the practical relevance of the current 
work. First of all, in many practical applications, it is preferable to use contrast agents 
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that do not emit photons after absorbing incoming excitation photons. In other words, 
the biochemical markers have quantum efficiency rj = 0. In this case, the second equation 
in ( 1 ) drops out of the transport system, and the terms involve 77 and all drop out 
from the datum (3). We are therefore back to the same mathematical problem as in a 
regular quantitative PAT process. The theory of the reconstruction in this case is covered in 
Theorem 3.3 of our results. Our results in this paper are in fact more general in the sense that 
we can deal with the general case of non-negligible quantum efficiency, that is 77 > 0. When 
77 > 0 , we have to take into account the impact of the emitted fluorescence photons in the 
reconstruction process. Neglecting this impact in the model would certainly introduce errors 
in the images reconstructed. The second issue we need to address is the difference between 
the work we have here and the theory on the same problem that have been developed in the 
diffusive regime [57] . It is generally believed that the radiative transport equation is a more 
accurate model than the diffusion equation to describe the propagation of NIR photons in 
biological tissues [9, 55], even though it is more complicated to theoretically analyze and 
numerically solve. Our analysis in this paper is useful when the diffusion approximation to 
the radiative transport equation breaks down, for instance in media of small volumes but 
large mean free paths. Optical imaging of small animals [33], for instance, is one of such 
biomedical applications for our work here. 

The rest of the paper is organized as follows. We first present in Section 2 some general 
properties of the inverse problem, especially the continuous dependence of the datum H 
on the unknown coefficients. We then consider in Section 3 the reconstruction of a single 
coefficient from a single internal data set. We derive some uniqueness and stability results 
on the reconstruction. In Section 4 we study the problem of reconstructing two coefficients 
simultaneously, mainly in linearized settings. We then present some numerical simulations 
based on synthetic data in Section 5 to validate the theory and the reconstruction algorithms 
we developed. Concluding remarks are offered in Section 6 . 

2 General Properties of the Inverse Problems 

We review in this section some general properties of the inverse problem of reconstructing 
77(x) and/or aa^xfi'^) in the transport system (1) from the datum H in (3). We denote 
by LP{X) (resp. the Lebesgue space of real-valued functions whose p-th power 

are Lebesgue integrable on X (resp. 12), and 'Hp{X) the space of L^(X) functions whose 
derivative in direction v is in U’{X), i.e. T-0p{X) = {/(x, v) : / G U’{X) and v • V/ G 
L^(A)}. We denote by L^’(r_) the space of functions that are traces of T-Lp{X) functions 
on r_ under the norm ||/||LP(r_) = /§rf-i |n(x) • w\\f\Pdwd'yf/P, d'f being the surface 

measure on dfl and = {v : v G s.t. — n(x) • v > 0}. It is well-known [2, 22] that 
both 'Hp(A) and iA(r_) are well-defined. To avoid confusion with 'Hp(X), we use 
to denote the usual Hilbert space of functions whose partial derivatives up to order 

k are all in Besides the assumptions in (A-i)-(A-ii), we assume further that: 

(A-iii) The domain is simply-connected with boundary dfl. The known optical coef¬ 
ficients satisfy 0 < Ci < cfa^xi, cra,m, crs,m,'^ < C 2 < -|-oo for some positive constants Ci 
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and C 2 . The unknown coefficients, {ri,aa,xf) belongs to the class 

A = {{r], cra,xf) : 0 < C 3 < r] < C 4 < 1 , 0 < C 5 < aa,xf < Ce < + 00 } (5) 


for some positive constants C3, C4, C5 and Cq. The scattering kernel 0 is symmetric, bounded 
and normalized in the sense that 


0(v,v') = 0(v', v), 0 < C7 < 0(v, v') < Cs < + 00 , Vv,v' e 

[ 0 (v,v')dv'= / 0 (v',v)dv'= 1 , Vve§‘^“\ 


d -1 


( 6 ) 


for some positive constants Cy and Cg. The illumination Pa;(x, v) is strictly positive such that 
0 < Cg < 5 - 3 ; (x,v) for some Cg. 

With the above settings, it is easy to see, following standard results in [2, 22], that the 
system ( 1 ) admits a unique solution in the following sense. 

Lemma 2.1. Let p G [l,oo] and assume that holds. Then for any given function 

gx{x,v) G L^’(r_), there exists a unique solution {ux,Um) & TLp{X) x TLp{X) to the couple 
transport system (1). Moreover, the following bound holds: 

||'fia:||LP(X) + ||Mm||LP(X) < c||pa; ||LP(r_) (7) 

with the constant c depending only on fl and the bounds for the eoefficients in assumption 

(A-iii). 


Proof. When the assumptions are satisfied, it follows directly from standard transport the¬ 
ory in [2, 22] that the first transport equation admits a unique solution Ux G 'Hp{X) such that 
||“a;||LP(x) < c||pa;||LP(r_). We then deduce, with the same argument that the second equa¬ 
tion admit a unique solution Um G Ttp^X) such that ||uto||lp(x) < ^\\Wa,xfKi{ux)\\LP{Q.) < 
c||mx||lp(x)- The bound in (7) then follows from selecting c = c(l -|- c). □ 

The above lemma ensures that the datum H in (3) is well-defined for any Pa;(x, v) G 
Z/’(r_) (p G [1, oo]) that satisfies the assumptions in (A-iii). Moreover H G LFiPl) following 
standard results in [22]. The next result shows that the datum H depends continuously on 
the unknown coefficients and is differentiable with respect to the coefficients in appropriate 
sense. 


Proposition 2.2. Letp G [1, oo] and assume that (A-iiiJ holds. Then for any given function 
Pa;(x, v) G Z/^’(r_), the datum H defined in (3), viewed as the map 


^a,xf\ 


iPi '^a,xf) 

L°°{n) X L°°(r2) 




LP{fl) 


( 8 ) 


is Frechet differentiable at any {r],(7axf) G x in the direction {Sr],5aaxf) G 

L^{n) X L°°{Q) that satisfy {j],<7a,xf) G A and {rj + Sr],aa^xf + ^<^a,xf) G A. The derivative 
is given by 

^ *^a,x/] rj)5Gaxf^Xj{Uxf)~\~0'^,j.Kj{Vx)~\~Crg^pjiKj{v.fff^ (9) 
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where {vx,Vm) G 'H^{X) x 'Hp(X) is the unique solution to 

V * XVx ~\~ e^t^x'^x ^s^xXq^Vx^ ^^a^xf^x^ X 

V ■ ^t,rn^m ^s,mXQ{Vjyi^ “1“ f]^a,xf Xj(^Vx^ (^^h(J(i^xf “1“ (^x)) X 

Vx{x,v) = 0, t;m(x,v) = 0 on r_ 

(10) 

where (Jt^x ^a,x 0~ ^s,x ond (Jt^m ~t~ ^s^m' 

Proof. Let f = 7] + Sr], aa,xf = (^a,xf + ^(^a,xf, and define X{r]aa,xf) = Wa,xf - Wa,xf- We 
denote by {ux, Um) the solution to (1) with the coefficients (77, ^a,xf)i and H the corresponding 
datum. It is straightforward to verify that (u^, u'^) = {ux — Ux, Um — Mm) solves the following 
system of transport equations 

V -Xu'^P at^xu'x = crs,xKeiu'J - 5aa,xfUx, in X 

V -Xu'^P = crs,mKe{u'^) Priaa,xfKi{u'^) P F{x.), in X (11) 

m(.(x, v) = 0 , m(„(x, v) = 0 on r_ 

with X(x) = X{r]aa,x)Xi{ux), and (n", u'ff) = [u'^ — Vx,u’^ — Vm) solves the following system 

V • Vm" + (Tt^xu'x = (Js^xKeiUx) ~ ^(^a,xfu'x, in X 

w -Xu'f^P at^mu'f^ = CFs,mKQ{u'ff) Pr]aa,xfKi{ul) PG{r^), in X (12) 

<(x,v) = o, u;;(x,v) = o on r_, 

with G'(x) = X{r]aa,xf)Ki{u'x) + Sr]Saa,xfKi{ux). 

With the assumptions on the coefficients and the illumination source Qx, we conclude 
that {ux,Um) e l-Lp{X) X 'Hp(X) and (fix,Mm) ^ Xfp{X) x 'Hp(X) [2, 22], This implies that 
F G and 

||-^||LP(f2) ||(^*^*^a,x T ^^^a,xf T ^V^^ci,xf^X j(^Ux^\\lp{Q,) 

— (Cl||(^fi||L°°(a) + ^2\\So'a^xf\\L^{a) + C3||(5r7||Loo(Q)||hcr(j_x/||L°°(r2))||fix||LP(X) 

— (ci||hr7||ioo(f2) + C211 hcTa^x/11 L°°(r2) + C3||hr7||/,oo(f2)||5cra,x/IU°°(a))||fi'x||LP(r_), (13) 

Following the same argument as in Lemma 2.1 we conclude that (11) admits a unique 
solution (m^,m(„) G Xp{X) x 'Hp(X) that satisfies 

||Mx||lp(x) < c||hcra,x/fix||LP(x) < c||hcra_x/IU°°(r2)||fix|| lp(x) < c||(5(Ja^x/IU°°(a)HwIUp(r_), (14) 
and 

||m(„||lp(x) < c{\\r)aa,xfXi{u'x)\\LP(a) + ||-^||LP(n)) < ^(||Mx||lp(x) + ll^llLP(a)) 

< (Cl||h?7|iL°°(n) + C2||hcra_x/||L°°(f2) + C3||hr7||ioo(Q)||(i(Ja^x/||L°°(Q))||fi'x||LP(r_)- (15) 
Therefore we have G G LP(fl) and the bound 

||G'||lp(Q) < \\{r]^(Ta,x + Sr](^a,xf + Sr]Scra,xf)Xi{u'^)\\LP(n) P \\Sr]S(Ta,xfXi{Ux)\\LPin) 

— (Ci||^?7lU°°(0) + C2||5cra^x/||L°°(Q) + C3||577 ||loo(q) ||(5cra^x/||L°°(f2)) ||m^||lp(Js:) 

+ ||^fi||L°°(a)||hcra,x/||L«>(n))||Mx||LP(X) 

< (c'Z 11 hr] 11^00(52) + C211 hcTa.x/11 L°°(f2) + C3 ||hr7||ioo(f2) ||hcra,x/|| l°°(Q)) HhcTa^x/IU°°(Q) Hw IUp(r_) • 

(16) 
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We then deduce, in the same manner as above, that (12) admits a unique solution 
that satisfies 

ll'lf^llLP(X) < ^\\S(Ja,xfU'J\LP(a) < ^\\S(7a,xf\\LP°(a)\\u'x\\LP(X) < Cc||(5cra^a,/Hf/x ||LP(r_), (W) 
and 


||“mlU»’W) ^ \\V(^a,xfKi{u")\\LP{n) + ||G'||LP(a) < c\\u"\\lp{x) + ||G'||LP(n) 

— (Cl||fi^^||L°°(n) + C2 IIhcTa x/11^00(^2) + C 3 ||<^^||l°°(Q) ||l°°(Q)) || j;^oo(q) ||LP(r_) • 

(18) 

The estimates (17) and (18) show that {u^, Um) is Frechet differentiable with respect to rj and 
(^a,xf as a map: L°°(fl) x L^{fl) i -4 L^(ll) x L^’(fl) {p e [1, oo]). Note that is independent 
of 77 , so its derivative with respect to p is zero, as can be seen from (17). 

The differentiability of H with respect to {p, <Ta,xf) then follows from the chain rule and 
the fact that is differentiable with respect to ( 77 , (Ta,xf)- Alternatively, it can also be seen 
easily from the bounds (14), (17), (18) and the following algebraic calculation: 


H[Pi H[p, 0'a,^xf] H [Pi CTa^xf^i^Vi ^'^a,xf') 

= E a2^^Ki{u”) + Oa^mKiiu'l^) + {5aa,xf - A{paa,x)Ki{u'^) - 5p5aa,xfKj{ux) 

This completes the proof. 


(19) 

□ 


We will study Born approximation, i.e. linearization, of the inverse problem of QfPAT 
in Section 4. The above result justifies the linearization process. To compute the partial 
derivative with respect to 77 (resp. (Ja,xf), denoted by H!^[p,aa,xf] (resp. H'^[p, aa,xf]), we 
simply set Saa^xf = 0 (resp. 5p = 0) in (9) and (10). It is straightforward to check that 


H^[V,<^a,xf]{Sv) 

^^a,xfKi (Ux') 


-5p + 


%,xfKj (Ux) 


A"/ {V-m) ; 


with Vm G 77 p(X) the unique solution to 


V • T 

'ym(x,v) 


T ^fi*^a,a;/A^/(^a:)) ill ^ 

0 , on r_. 


and 


e;(i - p)Ki{ux) 


= her, 


(TJ. 


a,xf 


(1 - p)Ki{Ux) 


Kl{Vx) 


(1 - p)Kj{ux) 




with {vx,Vm) G 'Hp(X) X 'Hp{X) the unique solution to 


V • XVx T CTtjx'^x ^s,a;A"©(u 3 ;) 5 (^a,xf'^x) in ^ 

^ ^T T V^ci,xfKli,^x^ T A"/(tij,), in X 

U 2 ,(x,v) = 0, 7;^(x,v) = 0 on r_. 


( 20 ) 


( 21 ) 


( 22 ) 


(23) 


The following result is a standard application of the averaging lemma [22, 28, 45]. It will 
be useful in Section 4. 
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Lemma 2.3. Assume that (A-iii ) holds. Let Qxi'x., v) G Z/°°(r_) be such that Ki{ux) > c > 0 


for some constant c. 
operator 


Then the rescaled linearized data 


viewed as the linear 


H(j \T ) ^a,a:/] {^'^a,xf ) _ dcTa^xf (1 'h)^^a,xf i(Vjn) 

EKi{ux) ■ L'^iQ) ^ L'^iQ) ^ ’ 

is Fredholm. The same is true for if the background coefficient (7a,xf > c > 0 

for some c. 

Proof. Let us denote by Sz {z G {x,m}) the solution operator of the transport equation 
with coefficients aa,z, (7s,z and vacnum boundary condition, i.e. Wz = Sz{f) with Wz the 
solution to: 

V •\/Wz + at,zWx - as,zKe{wz) = f, in X, Wz = 0 on r_. 

We can then write Kj{vx) and Kj{vm) in (24) respectively as 

-^/(w) Ax(^S(7a,xf') 1 and, Ki{yxn) •^ma:(^*^a,a;/) F ^m{dl^(7a,xf^ (2^) 

where the operators Aa,, and Kmx are defined as 

^xi.^^a,xf') — A"/(>S'a;('Ua;h(Tfl, a-j)), Axfi{b(7a,xf^ Xj(^Sm(,Kj(^Ux^S(7a,xf^^^ (26) 

A-mx(^*^a,x/) = ■kTl(^Sxn{j](7a,xf-kTl{^Sx{Ux5Ga,xf^^^- (2f) 

Following the averaging lemma [22, 28, 45] and the compact embedding of to Lf{Q), 

we conclude Kj : L‘^{X) —>■ is compact. Due to boundedness of Ux (and therefore 

Ki{ux)), T] and (7a,xf, both Sx and Sm are compact as operators from L^(D) to L‘^{X) with 
the assumptions on the coefficients in (A-i) [22, 45]. Hence, Aa;, A^, and A^a; are all 
compact operators on LA{Q). Therefore as an operator can be represented 

as (1 — rj)X + K, with K compact. Therefore it is Fredholm. The same argument works for 

^T7 I—I 

EKi{ux) • 

3 Reconstructing of a Single Coefficient 

In this section, we consider the reconstruction of one of the two coefficients of interests, 
assuming the other is known. We start with the reconstruction of the quantum efficiency. 

3.1 The reconstruction of r] 

Assume now that the fluorescence absorption coefficient (7a,xf is known and we are interested 
in reconstructing only t]. This is a linear inverse source problem. We can derive the following 
stability result on the reconstruction. 









Theorem 3.1. Letp G [1, cxd] and the source G L^iT-) he such that the transport solution 
Ux to (1) satisfies Kfiux) > c > 0 for any {r],aa,xf) ^ A. Let H and H be two data sets 
generated with coefficients {r],(7a,xf) o,nd {f],aa,xf) respectively. Then H = H a.e. implies 
g = fj a.e.. Moreover, the following stability estimate holds, 

c\\H — H\\lp{q) < II“ 'n)^a,xfKi{ux)\\LP(n) < C||i7 — H\\LP(n) (28) 

where the constants c and C depend on Ll and the coefficients aa^xi, <7a,m, <^s,x, and S. 

Proof. Let {ux,Um) and {ux,Um) be solutions to the coupled transport system (1) with 
coefficients {g, <Ta,xf) and {f}, (Ta,xf) respectively. We notice immediately that Ux = Ux- Define 
uJm = 'iJ'm ~ ^ra- We then verify that 

id')f^ (^g g^^a,xfLdj(^Ux^ 

This leads to the bound 

^ Cl||(^ V')^ci,xfLdi{^Ux')\\LP{Q,) T ^2II.^7(^m) ||lp(Q)• 

and the bound 

||(^ V^^a,xfLdj(^Ux^\\LP{Cl) — Cl(‘^)||.^ .^||lp(Q) T C2(*^a,m) ||.^7(^m) ||lp(Q)) 

We check also that solves the transport equation 

V • ^Wjyi T {p^a,m T (^rn) T (jl V^^a,xfLdl(Ux^ ^ 1^ W 

w^(x,v) = 0, on r_. 

It then follows from classical results in transport theory [2, 22] that this equation admits a 
unique solution Wm G TL^^X) that satisfies the following stability estimate 

||^m||LP(X) ^ *'3(D, (Ta jjj, (Tg jjj, .^) II ( 7 ^ g^^a,xfLdj(Ux^^L,P{Q)- (^3) 

The left bound in (28) then follows from (30) and (33). 

To derive the right bound in (28), we replace the last term in the transport equation (32) 
with (Ta,mKl{Wra) - (H - H)/E to get 

V • XWjyi T {p^a,m T ^a,mLdl(^Wxn} T ^s,mLdQ{Wrn) ) ln W (34) 

Wmixx) = 0, on r_. 


(29) 

(30) 

(31) 

(32) 


We define 0(x, v, v') = —-1-—0. It is straightforward to verify that 0 is 

symmetric and normalized in the sense of ( 6 ). We can then rewrite the above transport 
equation as 


V • VWjTj T ((Ta^m T ^s,m)^m 

M'm(x, V) 


ipa,m T ^s,m)L^ qI^Wyti) s ; in X 

0, on r_. 


(35) 


This is a transport equation for a conservative medium. Due to the fact that Q is bounded, 
classical results in transport theory (see for instance [22, Theorem 1 on page 337]) then 
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concludes that the equation admits a unique solution Wm G Moreover, we have the 

stability estimate 

||^m||LP(X) ^ ^a,mt ^s,mt ^^\\H -^||LP(f2) (26) 

The right bound in (28) then follows from (31) and (36). The uniqueness of the reconstruc¬ 
tion then follows from the fact that H = H implies = 0 from (35), which then implies 
r) = fj from (29). □ 

Note that the bound in (28) is weighted in the sense that it is on [r] — ff)Ki{ux) not 
directly on {r] — fj). This means that if Kiiu^) is too small, it is very hard to reconstruct 
accurately rj. 

The proof of the above stability result is constructive in the sense that it provides an 
explicit reconstruction procedure for the recovery of rj. We now summarize the procedure 
in the following algorithm. 

Reconstruction Algorithm I. 

51. Given cra,xf, solve the first transport equation in (1) with the boundary condition 
for 

52. Evaluate the function g(x) = <Ta,xKi{ux) — 

53. Solve the following transport equation for Um'- 

V • XUjji -\- (cTa^m ~\~ a,m T s,rri)K-\- (/(x), iu X f37) 

Um(x,v) = 0, on r_. ^ ^ 

54. Reconstruct rj as - - aa,xKi{ux) - <Ja,mKi{um))/{(Ja,xfKi{ux)). 

This is a direct reconstruction algorithm in the sense that it does not involve any iteration 
on the the unknown coefficient. The algorithm is very efficient since it requires solving the 
transport equation (37) only once. 

Remark 3.2. Thanks to the fact that the problem of reconstructing rj given Ca^xf is linear, 
we can easily verify that the same type of uniqueness and stability results in Theorem 3.1 
hold for the linearized problem of reconstructing rj defined in (21) and (20). Moreover, the 
above reconstruction algorithm works in exactly the same manner in the linearized setting. 

3.2 The reconstruction of cTa^xf 

We now assume that we know g and aim at reconstructing aa,xf- In this case, we can show 
the following result. 

Theorem 3.3. Let g^ G Z/^(r_) (p G [l,oo]J be such that the solution to the transport 
system (1) satisfies = Kjiux) > c > 0 for any coefficient pair {g,(7a,xf) G A. Let H and 
H be data sets generated with coefficient pairs ( 77 , (Ja,xf) and {rj, ^a,xf) respectively. Then 
H = H a.e. implies <Ta,xf = ^a,xf n-e.. Moreover, the following bound holds, 

C\\H — .^||lp(Q) < \\{(^a,xf — ^a,xf)Ki{Ux)\\LP{Cl) < C||i7 — H\\Lpjci), (38) 

with c and C depending on Ll, Ua^xi, ‘^s,x> 3 OL^d 5. 
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Proof. Let {ux,Um) and {ux,Um) be solutions to the coupled transport system (1) with 
coefficients {r],aa,xf) and {r],d'a,xf) respectively. Dehne Wx = Ux — Ux and Wm = Um — Um- 
Then we have 

^ ^ 

^ (^3^) T ^a,mKl{Wqxi) T (1 ^)(^a,x/ ^a,xf^Kl{Px) • (3®) 

This leads to the bound 

||.^ .^||LP(f2) ^ II.^7(^x) ||LP(f2) T ^2IIIIffi ^3 II *^a,x/).^7(^x) ||LP(f2); (^0) 

and the bound 


||(*^a,x/ ^a,xf^^l{Px')^LP{Q,) — ^ll|-^ .^||LP(f2) T *-2 ll"^.^ ll^'^’l^^) ffi ^3 ll"^.^ l|i'*’(i^)' 

We now verify that {wx,Wm) solves the following transport system: 

V • WWx T ^t,x'^x ^s,x^q{P^x) iPa,xf ^a,x/)^X) ffi ^ 

V • ^Wxn T CFt,m'^m ^s,m^Qip^m) T V^a,xf^li,^x^ T V(.^a,xf ^a,x/)-^7(^x)) ^ ^ 

Wx{:>^,v) = 0, tt;m(x,v) = 0, on r_ 

(42) 

where at^x = <^a,xi + ^a,xf + <^s,m- We then deduce, following similar procedure as in the proof 
of Proposition 2.2, that 


^x||lp(X) T ||^m||LP(X) — *-4ll(*^a,x/ *^a,x/).^7 (^x) || LP(f2) • 


(43) 


The left bound in (38) then follows from (40) and (43). 

To derive the right bound in (38), we use (39) to eliminate the quantity (Ja,xf — ^a,xf in 
the transport system (42) to obtain: 


V ■'\/Wx + at,xWx 
V • V'iCjjj T 

Wx{x,v) = 0 , 


as,xKe{wx) + a'^^Ki{wx) + cr; 


/ 

s^xm 


as,mKe{w^) - a'^^^Ki{wJ) - a'^ „^^Ki{wx) + 
u;^(x,v) = 0 , 


K^jl'ijn i — —— H)ux — 

(H-H)n 
v) ’ 


in X 
in X 

on r_ 

(44) 


where crl ^ = 


r htcc 


s.x (l-r])Kj(ux)^ ^s,xm {l-ri)Ki{ux)'> ^s,m 


(Pa,mUx 


(J. 


_ 'n^a,'n 

l-T] 


and a' 


_ 'n^a,x 

l-T] 


To write the 


system in standard form, we perform the change of variable Wx —Wx- We then have 


v-Vwx + at,xWx + af^^Kj{wm) = (7s,xKe{wx) + a'^^^Kj{wx) + ^ 

V • ^Wxn T T ^s,m^ 1 T ^s,mx^I^'^x) T ^S(l—?;) ’ i^ W 

Wx{yi,^) = 0 , M;m(x,v) = 0 , on r_ 

(45) 

With the assumption on the coefficients cr^ j., cr's^xm^ ^s,mx positive. We 

check also, after using the assumption Ux = Kjiux), that Ai = dt^x + <^'s,xm ~ W,x — <^'s,x — 
da,x + (cra,m-d2,x)/[(l-^)] ^nd A 2 = CTj,™ + C^'“ W,m “ C^' ^mx = (c^a.m “ Wa,xi)/(1 “ ??) • The 
conditions in Theorem 3.3 ensure that Ai, A 2 > c' > 0 for some d. We can therefore combine 
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the techniques in [29, 63, 64, 71], see detailed analysis in [53], to show that system (45) admits 
a unique solution {wx,Wm) that satisfies 

||Wa,||LP(x) + ||Wm||LP(X) < — -^||LP(a)- (46) 

We can now combine (41) and (46) to obtain the right bound in (38). The uniqueness 
result follows from the fact that (45) admits only the trivial solution [wx, Wm) = (0,0) when 
H = H. □ 

Linearized Case. Unlike in the case of reconstructing rj, the above proof is not construc¬ 
tive since the unknown coefficient (Ja,xf shows up in the transport system (45). Therefore, 
the proof does not provide directly a reconstruction algorithm. For numerical reconstruc¬ 
tions in this nonlinear setting, we use the optimization-based algorithm in Section 4.4. If 
we consider the same problem in linearized setting, we can indeed derive an explicit recon¬ 
struction procedure. To do that, we replace the Saa,xf in (23) with its expression given in 
the linearized datum ( 22 ) to get the following system: 

V ■ '’^Vx 0's,xmKl{Vm) = 0 's^xKq{Vx) + <Jg x^l{'^x) ~ {l-rf)'E.Ki(ux)' ^ 

^ ■^Vm + <yt,mVm +Cr'^,mKi{Vm) = as,mKQ{Vm) + i{Vx) + mX 

Vx{yi,v)=0, v^(x,v)^0, on r_ 

(47) 

where we have performed the change of variable Vx —>■ —Vx, and the coefficient = 

7 ? 

) ■> while the coefficients (j'gxm^ ^s,m-> ^'s,mx defined as in (44). This system 
does not contain the unknown coefficient boa^xf- It can be solved for {vx,Vm)- We can then 
reconstruct 5(7a,xf following ( 22 ). The reconstruction procedure can be summarized into the 
following reconstruction algorithm. 

Reconstruction Algorithm II. 

51. Given the background coefficient (Ta,xfi solve the first transport equation in ( 1 ) with 
the boundary condition Qx for Ux (and therefore Ki(ux))', 

5 2 . Evaluate the coefficients a's^xm^ K,m 

53. Solve the transport system (47) for {vx,Vm) and perform the transform {—Vx,Vm) —t 
(W) ^m)) 

H' 

54. Reconstruct 5(Ta,xf as [^ - a^^^Ki{vx) - cr<j,^A/(n^)] / [(1 - ri)Ki{ux)]. 

Following the control theory for transport equations developed in [1, 3, 38], we can show, 
under reasonable assumptions, the existence of sources Qx such that Ux = Kj{ux) holds for 
each pair ( 77 , cra,xf) £ Such sources, however, might be complicated, for instance we 
might need to solve a control problem, to construct in practical applications. The usefulness 
of Reconstruction Algorithm II is therefore limited by this fact. Note that in applications 
where the medium is scattering-free, see for instance discussions in [23, 44], this algorithm 
is indeed very useful since there are many ways to construct illuminations sources to have 

Ux Kj(Ux^. 
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4 Simultaneous Reconstruction of Two Coefficients 

We now consider the problem of simultaneous reconstruction of the quantum efficiency and 
the fluorescence absorption coeflicient. We start with the linearized case. 

4.1 Linearization around {r],(Ta,xf) — (0,0) 

We first consider the special case where both coefficients are small. In this case the product of 
the coefficient is small so that generation of fluorescence is very small and can be neglected. 
Therefore, the system involves only the light at the excitation wavelength. The QfPAT 
problem reduces to the usual quantitative PAT problem. To be precise, we linearize the 
problem around the background {r],aa,xf) = (0,0). Then the second transport equation 
in (10) has the solution = 0. Therefore, the datum (9) simplifies to 

[O, O] (hr^, (5 cT(j j-y) hcT(j A"/(Uj,) A ^a,xiKj{Vx)i (18) 


and the first transport equation in system (10) simpiifles to 

V • WVx {(^a,xi T *^s,a:)W *^s,a:A^0(w) 111 

?;a;(x,v) = 0, on r_. 

We observe that 5r] does not appear in the datum (48) or the equation (49). Therefore, it 
can not be reconstructed in this setting. We can show the following result. 

Proposition 4.1. Let Ux be the solution to the first transport equation m (1) with cra,xf = 0. 
Let Qx G Li’(r_) (p e [l,oo]) be such that Ux = Kfiux) > c > 0. Denote by iif'[0,0] and 
H'[0, 0] the perturbed data sets in the form of (48), generated with perturbed coefficients (Sr], 
b<Ja,xf) CLnd (5r], 5aa,xf) respectively. Then iif'[0,0] = i7'[0,0] a.e. implies 5aa,xf = Saa,xf 
a.e.. In addition, we have, 

c||i/'[0,0]-^'[0,0]|Up(o) < \\iSaa,xf-^xj)Kiiux)\\Lnn) < C||//'[0,0]-^'[0,0]|Up(o), (50) 
with c and C constants that depend on Ll, S, (Ta,xi o,i^d ag^x- 
Proof. The datum (48) implies directly that 

II [0,0] H [O, O] ||lp(Q) ^ Cl II {SCaxf b(7 (ixf)KI {Ux) ||lp(Q) T *-21| W W ||lp(X) , (^1) 

and 

\\{S(ya,xf — baa,xf)Ki{ux)\\LP(n) < c'^||-^^[0,0] — H '[ 0 , 0] ||lp(q) + — hx\\LP(x), (52) 

with the constants depend on fl, aa^xi and S. 

With the assumptions in the theorem, we deduce from the transport equation (49) that 

||W W||lp(Q) — ^3II (^ll'ajX/ bo'a^xf)^fi^x')\\LP{Cl)- (^8) 
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The left bound in (50) then follows from (51) and (53). To get the right bound in (50), we 
use the datum (48), and the assumption that Ux = Ki{ux), to rewrite (49) as 

^ ■ '^'^x T iS^a,xi T O'g x^Vx (^s,x^o{,^x') T ^a,xi^li,^x) 5 ; 1^^ ^ f54') 

Vx{x,v) = 0 , on r_. 

This is a conservative transport equation that admits a unique solution with the stability 
result: 

\\vx - VxWl^u) < c'3||i/'[0,0] - //'[0,0 ]|Up(o), (55) 

where Cg depends on 0, cTa.cco (^s,x and E. The right bound in (50) then follows from (52) 
and (55). ’ ’ □ 

The above proof is again constructive when a gx that satishes the assumption in the 
theorem is available to us, in the sense that we only need to solve (54) for Vx and then 
compute Saa,xf = {H'[0,0]/E - aa,xiKi{vx))/Ki{ux). 


4.2 Linearization around a general background 

We now consider the linearization around a general background (rj ^ 0,aa^xf ^ 0). We 
study the case where we have J >2 data sets, 1 < j < J: 


H'[ri,aa^xf]{Sv,Saa^xf) 


= i-Sricra,xf + (1 - v)Scra,xf) + 


cr'‘ 


Kjiui) 


Kj{vi) + 


Kiiui) 


Ki{vl) 


(56) 


where ul is the solution to the hrst transport equation in (1) with background coefficient 
(7a,xf and illumination source gl, while (u^,u^) is the solution to the coupled system (10). 

To study the linear inverse problem dehned in (56), we introduce two new variables 
C = Sgaa,xf + V^<7a,xf and ^ = Saa,xf- If is straightforward to verify that (C,0 uniquely 
determines (Srj, Saa,xf) when g ^ 0 and (7a,xf ^ 0. We can then collect the J data sets to 
have the following linear system for the unknown coefficient pair (C,0- 


n 


/ -X + X - nj \ 


= z, with, n = 


\ -X + x-W 


i J 




and z = 


=.Ki{ul) 




(57) 


with n 


j _ <Ta,m 

C “ Kiiui) 


and 


Tj 

Kiiui) ^ ^ Kiiui) 


Here A^, A^^. and A^ are defined 


as in (26) and (27) with Ux replaced by From Lemma 2.3 we know that H;^ and H^ 
(1 < j < J) are compact operators on 

From the discussion in the previous sections, we know that X — H^ and X — n| are all 
invertible for well-selected illumination sources g^, 1 < j < J. However, that does not 
guarantee the invertibility of the linear system (57). For the case of J = 2, the invertibility 
of the system (57) is equivalent to the invertibility of (X—n^)“^(X—n|) —(X—n^)“^(X—H^). 
Therefore, we need to choose illumination sources gl and g^ such that (X — n^)“^(X — n|) — 
(X- nj)-i(x- Hj ) is invertible; see next section for some discussions on the regularized 
version of this problem. 
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4.3 A partially linearized model 


We now briefly discuss a very popular simplification of the mathematical model in the flu¬ 
orescence optical tomography literature. This simplification assumes that the fluorescence 
absorption coefiicient <7a,xf is small compared to the background tissue absorption coeffi¬ 
cient (Ta^xi- Therefore, it can be dropped from the first equation in the model (1); see for 
instance [46]. In other words, the model, for source gl {1 < j < J), now reads, 


V • Vul -s {aa,xi + <ys,x)ul 

3,770 (li^), 

in 

w 


V ■ -S {aa,m + crs,m)ui. 

= + gaa,xf Kiiui), 

in 

w 

(58) 

<(x,v) =9l, 

<(x,v) = 0 

on 

r_. 



The data, for source gl {1 < j < J), now simplify to, 


H, = 




7 / 1 \ ^ • 

^a,xi (1 V)^a,xf T 


Ki{ui) 


Kiiui) 




(59) 


The inverse problem of reconstructing rj and (Ta,xf from datum (59) is a nonlinear problem 
despite the fact that a partial linearization has been performed on the transport model. 
However, if we define ((^ = (1 — ri)<Ja,xf and ^ = Ua^xf, then the inverse problem is bilinear 
with respect to {C,0- Precisely, we can write the datum as, 


n, = (- nj(C) + nj(«). 1 < i < J 


(60) 


with H;^ = ^ defined the same way as before and being compact on This can 

again be written into the form of linear system (57) with the coefficient matrix and source 
vector respectively 



/ X - nj nj \ 


( Hi \ 

n = 


, and, z = 





\Hj ) 


(61) 


Regularized Inversion with J = 2. In the case that two data sets are available, we can 
solve the inverse problems in this section and Section. 4.2 in regularized form. To do that, 
we observe that if we define 

n„ = n+(° a >0 (62) 

then Ha is a Fredholm operator on x for the 11 defined in both (57) and (61). 

To be precise, 11^ are respectively. 


n 


a 


-i + ul 
-x + n2 


I 




aX + X-Ul 


-x + ul 
n^-nj 


X 




aX + Yi\- n| 


-x + nj 


X 

aX + 


+ 


n^-n^ 



(63) 
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and 


(64) 


^(T-ul ul 0 \ / 0 nn 

" \i-ul ai + ui J V ^ / V -nj nj y 

where ~ is used to denote the elementary operation of subtracting the first row from the 
second row. For any fixed a > 0, let us denote by A/'(na) the null space of matrix operator 
IIq , then the following result follows immediately from classical stability theory of Fredholm 
operators [35]. 


Proposition 4.2. Let z and z be two perturbed data sets defined as in (57) or (61). Let 

C ^ f C 


(C)0* solution to n, 

a > 0. Then we have 




= z and Ilr 




= z respectively for some 


c||z - z||(i 2 (n ))2 < ||(c,0 - (C,0ll(L2(Q))2/A/-(n„) < C||z - z||(i 2 (n)) 2 . (65) 


for some constants c and C. 

In the numerical computation, to solve (57) or (61) directly, we have to construct the 
operator 11 explicitly. This is hard to do in practice since it essentially requires the analytical 
form of the Green’s function for the transport equation at the emission wavelength. We do 
not have access to this Green’s function. Instead, solve the linear problem with a classical 
method of Landweber iteration [36] that we summarize in the following algorithm. 

Reconstruction Algorithm III. 

51. Take initial guess (Co)^o); 

5 2 . Iteratively update the unknown through the iteration: 


f Ck+l 

\ C/c+l 


(I - rn*n) 


Ck 

Ck 


+ rn*z, k > 0. 


( 66 ) 


Stop the iteration when desired convergence criteria are satisfied. 


Here r is a positive algorithmic parameter that we select by trial and error, 
operator 11 * is formed by transposing 11 and replacing H^ and H^ with H^* = 


Kr 


Kiiui) 


and nf = uis; 


oKt 


Kiiui) 


+ KSl oKjO oKjO 


Kiiui) 


The adjoint 

Ki(4)S-„ o 

respectively. 


Here S'* is the adjoint of Sz {z G {x, m}) that is defined as the solution operator of the 
adjoint transport equation with coefficients aa,z, <^s,z and vacuum boundary condition, i.e. 
Wz = S*{f) with Wz the solution to: 


-v-Vwz + {(Ta,z + crs,z)w^-(7s,zKe{Wz) = f, in A, = 0 on r+. 


Therefore, at iteration k of the Landweber algorithm, we solve J forward transport systems 
and then J adjoint transport systems to apply the operator 11*11 to the vector {Ck^^kY- 
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4.4 Iterative reconstruction for the nonlinear case 

For the simultaneous reconstruction of r} and (7a,xf in the general nonlinear case, we do not 
have any theoretical results on uniqueness and stability currently. Nor do we have more 
explicit reconstruction methods. We rely mostly on general computational optimization 
techniques to solve the inverse problem. More precisely, we search for solutions to the 
inverse problem by minimizing the objective functional: 

(7a,xf) = 2 ^ dx + ^R{ri, aa,xf) (67) 

i=i 

where the regularization functional is taken as R(ri, (7a,xf) = |(|| Vr]||^^ 2 (Q)]d+|| Vcra,a;/||^^ 2 (Q)]d)- 
Following the result in Proposition 2.2 and the chain rule, we can obtain the following 
result straightforwardly. 

Corollary 4.3. The functional ^{'q^(7a,xf), viewed as the map: 4 : Wf{VL) x i-?- E+ 

is Frechet differentiable at any {ri,aa,xf) G ^ F\ A. The partial derivatives in 

the direction 5r] (such that {r] + Sr],aa,xf) ^ A) and the direction Saa,xf (such that {'r],aa,xf + 
^vra,xf) & A) are given respectively as 


KVl.(^a,xf]{5vi) = 


$ 


- Sr]aa,xfKi{ul.) + aa,mKi{wl^)] + /3VSr] ■ Vr^jd^S) 

a[V^ ((a,xf]iS(7a,xf) ~ I ^ ] ^j^\^(7a,xf(^ ~ (^a,x^li^x) W,rn-77j('^m)] hx(69) 

da ..-I 

(70) 


i=i 


+ 13 VSaa,xf ■ Vaa,xfdy-, 
Ja 


where the residual Zj = 'E[af^Ki(ui.) + aa,mRi(uf()] ~ Hj, Wm is the unique solution to (21), 
and {vx,Vm) is the unique solution to (23). 


We can therefore employ gradient-based minimization techniques to minimize the func¬ 
tional (67). Here we use the limited memory version of the BFGS quasi-Newton method 
that we implemented in [54]. This method requires only the gradients of the objective func¬ 
tional which we derived in Corollary 4.3. To simplify the computation of these gradients 
numerically, we apply the adjoint state technique. We denote by (g^, qlf) the unique solution 
to the following adjoint transport system: 


-V • Vqi + at,xqi 

-V • Vg 4 + (^t,mqL 

?^(x,v) = 0 , 


(^s,xKe{qi) + ^(^2,x^j + Wa,xfKi{qff), in X 
as,mKe{q(a) + '^vra,mZj, in X 

<(x,v) = 0 on F+ 


(71) 
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It is then straightforward to show that 


K^V.(^a,xf]{Sv) 


r ^ 

= / {E*'<'«-f'f'^'K)[-S5,+A-,(9^)]+/JVi,.V.,}dx,(72) 

Jn 

= / ^^ci,xfKi{ul.) [S(l - r))zj + r]Ki{q^) - Ki{q^)] dx 

Jn 

+ 13 V5aa,xf • S/Ga,xfdyi. (73) 

Jq 


Therefore, to compute gradients of the $ at {ri,aa,xf)i we only need to solve a set of J 
forward transport systems (1) and a set of J adjoint transport systems (71). We can then 
evaluate the gradients in any given direction {5r],Saa^xf) according to (72) and (73). 

It is obvious that this optimization-based nonlinear reconstruction method can be used 
also to reconstruct a single coefficient. To only reconstruct rj, we only need to set the gradient 
with respect to Ga^xf fo zero and vice versa. 


5 Numerical Experiments 

We now present some numerical reconstructions using synthetic interior data. We restrict 
ourselves to two-dimensional settings only to simplify the computation. 

The spatial domain of the reconstruction is the square 12 = (—1,1) x (—1, 1). All the 
transport equations in 12 x are discretized angularly with the discrete ordinate method 
and spatially with a first-order finite element method on triangular meshes. In all the 
simulations in this section, reconstructions are performed on a finite element mesh consisting 
of about 2000 triangles and a discrete ordinate set with 64 directions. For the absorption 
and scattering coefficients that are known, we take 

Ga,xi = Ga,m = g’'^ {2-{l2x\ + l2y\ mod 2)), (74) 

Gs,x = Gs,m^Gl{l + {[2x] + [_2y] mod 2)), (75) 

where [-J represents the floor operation, and g^ are respectively the base level absorption 
and scattering coefficients. In all the cases below, we set cr^ = 0.1. The value of g^ varies from 
case to case and will be given below; see Fig. 1 (i) and (ii) for plots of the two coefficients. 
The scattering kernel 0 is set to be the Henyey-Greenstein phase function [9, 32, 70] which 
depends only on the product v • v'. 

To generate synthetic data for the nonlinear inversions, we solve the transport system (1) 
with true quantum efficiency r] and fluorescent absorption coefficient Ga^xf and compute 
H according to (3). To generate synthetic data for linearized inversions, for instance in 
Experiment 3 below, we use directly the linearized data models, for instance (56), with the 
true coefficient perturbations. This way, we can exclude the linearization error from the 
data used in the inversion. We do this since our main aim is to test the performance of 
the reconstruction algorithms, not to check the accuracy of the linearizations. To mimic 
noisy measurements, we add additive random noise to the synthetic data by multiplying 
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each datum point by (1 + 7 x 10 “^noriiirnd) with normrnd a standard Gaussian random 
variable and 7 a number representing the noise level in percentage. When 7 = 0, we say the 
data are noise-free. 

To measure the quality of the reconstruction, we use the relative error. This error is 
defined as the norm of the difference between the reconstructed coefficient and the true 
coefficient, divided by the norm of the true coefficient and then multiplied by 100 . 

We performed numerical simulations on the reconstructions of many different coefficients 
pairs ( 77 , o'a,xf)- The qualities of the the reconstructions are very similar. To avoid repetition, 
we will present only reconstructions for a typical coefficient pair we show in (iii)-(iv) of Fig. 1. 








I 



Figure 1: From left to right are: (i) the absorption coefficient aa^^i = <^a,m defined in (74) 
with = 0.1, (ii) the scattering coefficient ag^x = <^s,m defined in (75) with = 2.0, (iii) 
the true quantum efficiency 77 to be reconstructed in the numerical experiments, and (iv) 
the true fluorescence absorption coefficient (Ja,xf to be reconstructed. 


Experiment 1. In the first set of numerical studies, we consider the reconstruction of the 
quantum efficiency 77 assuming that the fluorescent absorption coefficient (Ta,xf is known. 
We use the Reconstruction Algorithm I presented in Section 3.1. We first perform numeri¬ 
cal experiments in isotropic medium with two different strengths of scattering coefficients. 
We show in Fig. 2 the reconstructions of 77 under base scattering — 1.0. Shown from 
left to right are respectively the 77 reconstructed using data with noise level 7 = 0, 2, 5 
and 10 respectively. The relative errors in the reconstructions are respectively 0.01%, 
14.24%, 35.59% and 71.18%. We repeat the simulations for a medium with stronger (but 
still isotropic) scattering (cr^ = 9.0). The results are shown in Fig. 3. The relative errors 
in this case are 1.04%, 14.84%, 37.02% and 74.02% respectively. If we compare the results 
in Fig. 2 and those in Fig. 3, we see that the quality of the reconstructions are almost inde¬ 
pendent of the scattering strength. This is what we observed in our numerical experiments 
in other cases as well. 

Experiment 2. In the second set of numerical studies, we consider the reconstruction of 
the fluorescent absorption coefficient (Ta,xf assuming that the quantum efficiency 77 is known. 
Currently, we do not have a well-established method to construct illuminations sources 
such that the condition Ux = Ki{ux) is satisfied for the transport solution, besides in non¬ 
scattering media. We therefore can not use directly the Reconstruction Algorithm II as we 
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Figure 2 : The quantum efficiency 77 reconstructed with different types of data. The noise 
levels in the data used for the reconstructions, from left to right are 7 = 0, 2, 5 and 10 
respectively. The base scattering strength is = 1.0. 
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Figure 3: Same as in Fig. 2 but with base scattering strength = 9.0. 


commented before. Instead, we use the nonlinear reconstruction algorithm in Section 4.4. 
We show in Fig. 4 the reconstructions of (7a,xf in an isotropic medium with base scattering 
strength = 1.0. Shown from left to right are respectively the reconstructions using data 
with noise levels 7 = 0 , 2 , 5 and 10. The relative errors in the four reconstructions are 
0.01%, 6.42%,16.06% and 32.12% respectively. In Fig. 5, we show the same reconstructions 
in an anisotropic scattering medium with base scattering strength = 9.0 and anisotropic 
factor 0.9. The relative errors are 0.02%,6.70%,16.74% and 33.42%, respectively. We 
again observed that the reconstructions are of good quality with data contains reasonably 
low level of random noise. 



Il 



Figure 4: The fluorescence absorption coefficient (7a,xf reconstructed with different types of 
data. The noise level in the data used for the reconstructions, from left to right are: 7 = 0 
(noise-free), 7 = 2, 7 = 5, and 7 = 10. The base scattering strength is = 1.0. 
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Figure 5: Same as in Fig. 4 but in a medium of anisotropic scattering with base scattering 
strength = 9.0 and anisotropic factor 0.9. 


Experiment 3. In the third set of numerical simulations, we study the simultaneous 
reconstruction of the coefficients rj and (Ta,xf in the linearized setting described in Section 4.3 
using the Reconstruction Algorithm III. The synthetic perturbed data are generated using 
directly the linearized model (56), not the original nonlinear model. Our aim here is to 
test the stability of the reconstruction, not the accuracy of the linearization. We use data 
sets collected from four angularly-resolved illuminations supported respectively on the four 
sides of the boundaries of the domain, pointing toward the interior of the domain. The 
background scattering strength is = 1.0 and the anisotropic factor is 0.5. We linearize 
the problem around the background coefficients: 

The reconstructions, after adding back the background, are shown in Fig. 6 . The refative 
error in the reconstructions using data with noise fevel 7 = 0, 7 = 2, 7 = 5 and 7 = 10 
are respectivefy (0.00%, 0.00%), (14.65%, 7.45%), (37.28%, 18.77%) and (75.80%, 39.04%) 
respectively. In all reconstructions, we applied the Tikhonov regularization with a small 
regularization strength that we select by trial and errors. We hope to develop more system¬ 
atical strategy on regularization in the future. 

Experiment 4. The last set of numerical simulations are devoted to the simultaneous 
reconstructions of the coefficient pair ( 77 , (Ta,xf) in the fufly nonfinear setting. We use the 
optimization-based reconstruction algorithm developed in Section 4.4. Besides the fact that 
the synthetic data are now generated with the full transport model ( 1 ), not the linearized 
model (56), the setup (for instance the background coefficients and anisotropic factor etc) 
is the same as that in Experiment 3. We performed reconstructions with data containing 
various noise levels. When the noise level is too high, we have difficulties to find reasonable 
initial guesses to make the algorithm converge. We show in Fig. 7 reconstructions with 
data containing a small amount of noise, 7 = 0 , 1 and 2 respectively, with the initial guess 
(^°Wa,a;/) being the average of the true coefficients inside the domain. The relative error in 
the reconstructions are respectively (16.40%, 8.32%), (18.26%, 9.17%) and (23.26%, 19.30%) 
respectively. We again impose weak Tikhonov regularizations in all the reconstructions with 
the regularization strengths selected by trial and error. Tuning various parameters in the 
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Figure 6: Simultaneous reconstructions of the coefficient pair (rj, <Ta,xf) in the linearized 
setting with different types of data. The noise level in the data used for the reconstructions 
are (from left to right): 7 = 0, 2, 5 and 10 respectively. The base scattering strength is 
al = 1 . 0 . 

algorithm could potentially improve the reconstructions results, but we did not pursue in 
that direction. 


6 Concluding Remarks 

We studied in this work a few inverse problems in quantitative fluorescence photoacoustic 
tomography in the radiative transport regime. We derived some uniqueness and stability 
results on the reconstruction of the fluorescence absorption coefficient and the quantum 
efficiency of the medium. In some cases, we were also able to develop efficient numerical 
reconstruction algorithms. These results complement the results in [57] for the QfPAT 
problem in the diffusive regime. We showed numerical simulations based on synthetic data 
to support the mathematical analysis and demonstrate the performance of some of the 
reconstruction algorithms. 

One important application of the results in this paper is in X-ray modulated fluorescence 
tomography (or X-ray luminescence tomography (XLT)) [61]. In XLT, X-rays, instead of 
NIR photons, are used to excite the molecular markers. The X-ray density and the 
generated NIR photon densities Um solve the coupled transport system (1) with the scattering 
term Kq{ux) = 0 since X-rays travel in straight fines without being scattered. The theory 
and reconstruction methods we developed in this work remain valid in that case. In other 
words, we can recover stably the fluorescence absorption coefficient using data collected 
from one X-ray illumination. This would provide a useful alternative to the reconstruction 
method for XLT in [61]. 

Even though the QfPAT problem has been analyzed in detail in [57] in the diffusive 
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Figure 7: Simultaneous reconstruction of the coefficient pair (rj, (Ta,xf) in the nonlinear setting 
with different types of data. The noise level in the data used for the reconstructions, from 
left to right, are respectively 7 = 0 , 1 and 2 . 


regime, the developments in this work are still useful in many settings. One well-known 
example is the application in optical imaging of small animals [33] where the diffusion model 
is not sufficiently accurate to describe the propagation of NIR photons inside the animals. 

Our main research focus in near future is to analyze the uniqueness and stability proper¬ 
ties of the simultaneous reconstruction problem, i.e. the problem of reconstructing the pair 
{v^^a,xf), in the fully nonlinear setting. This is an unsolved problem even in the diffusive 
regime [57], although numerical simulations we have so far suggested that uniqueness and 
stability both hold, at least in the regime where both coefficients are sufficiently large. 


Acknowledgments 

We would like to thank the anonymous referees whose comments help us improve signifi¬ 
cantly the quality of this paper. This work is partially supported by the National Science 
Foundation through grant DMS-1321018, and the University of Texas at Austin through a 
Moncrief Grand Challenge Faculty Award. 


References 

[ 1 ] S. Acosta, Time reversal for radiative transport with applications to inverse and con¬ 
trol problems^ Inverse Problems, 29 (2013). 085014. 

[ 2 ] V. Agoshkov, Boundary Value Problems for the Transport Equations, Birkhauser, 
Boston, 1998. 


23 





















[3] V. I. Agoshkov and C. Bardos, Optimal control approach in inverse radiative 
transfer problems: The problem on boundary function, ESAIM: Control, Optimisation 
and Calculus of Variations, 5 (2000), pp. 259-278. 

[4] M. Agranovsky, P. Kughment, and L. Kunyansky, On reconstruction formu¬ 
las and algorithms for the thermoacoustic tomography, in Photoacoustic Imaging and 
Spectroscopy, L. V. Wang, ed., CRC Press, 2009, pp. 89-101. 

[5] D. Alvarez, P. Medina, and M. Mosgoso, Fluorescence lifetime imaging from 
time resolved measurements using a shape-based approach, Optics Express, 17 (2009), 
pp. 8843-8855. 

[6] H. Ammari, E. Bossy, V. Jugnon, and H. Kang, Mathematical modelling in 
photo-acoustic imaging of small absorbers, SIAM Rev., 52 (2010), pp. 677-695. 

[7] H. Ammari, E. Bretin, V. Jugnon, and A. Wahab, Photo-acoustic imaging 
for attenuating acoustic media, in Mathematical Modeling in Biomedical Imaging II, 
H. Ammari, ed., vol. 2035 of Lecture Notes in Mathematics, Springer-Verlag, 2012, 
pp. 53-80. 

[8] H. Ammari, J. Garnier, and L. Giovangigli, Mathematical modeling of fluores¬ 
cence diffuse optical imaging of cell membrane potential changes, Quarterly of Applied 
Mathematics, (2012). 

[9] S. R. Arridge, Optical tomography in medical imaging, Inverse Probh, 15 (1999), 
pp. R41-R93. 

[10] G. Bal, Hybrid inverse problems and internal functionals, in Inside Out: Inverse Prob¬ 
lems and Applications, G. Uhlmann, ed., vol. 60 of Mathematical Sciences Research 
Institute Publications, Cambridge University Press, 2012, pp. 325-368. 

[11] G. Bal, A. Jollivet, and V. Jugnon, Inverse transport theory of photoacoustics, 
Inverse Problems, 26 (2010). 025011. 

[12] G. Bal and K. Ren, Multi-source quantitative PAT in diffusive regime, Inverse Prob¬ 
lems, 27 (2011). 075003. 

[13] -, On multi-spectral quantitative photoacoustic tomography in diffusive regime, In¬ 

verse Problems, 28 (2012). 025010. 

[14] G. Bal and G. Uhlmann, Inverse diffusion theory of photoacousties. Inverse Prob¬ 
lems, 26 (2010). 085010. 

[15] P. Beard, Biomedieal photoacoustie imaging. Interface Focus, 1 (2011), pp. 602-631. 

[16] P. Burgholzer, H. Grun, and A. Sonnleitner, Photoaeoustie tomography: 
Sounding out fluorescent proteins, Nat. Photon., 3 (2009), p. 378379. 


24 



[17] P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf, Exact and ap¬ 
proximative imaging methods for photoacoustic tomography using an arbitrary detection 
surface, Phys. Rev. E, 75 (2007). 046706. 

[18] B. T. Cox, S. R. Arridge, and P. C. Beard, Photoacoustic tomography with a 
limited-aperture planar sensor and a reverberant cavity, Inverse Problems, 23 (2007), 
pp. S95-S112. 

[19] B. T. Cox, S. R. Arridge, K. P. Kostli, and P. C. Beard, Two-dimensional 
quantitative photoacoustic image reconstruction of absorption distributions in scattering 
media by use of a simple iterative method, Applied Optics, 45 (2006), pp. 1866-1875. 

[20] B. T. Cox, J. G. Laufer, and P. C. Beard, The challenges for quantitative 
photoacoustic imaging, Proc. of SPIE, 7177 (2009). 717713. 

[21] B. T. Cox, T. Tarvainen, and S. R. Arridge, Multiple illumination quantitative 
photoacoustic tomography using transport and diffusion models, in Tomography and 
Inverse Transport Theory, G. Bal, D. Einch, P. Kuchment, J. Schotland, P. Stefanov, 
and G. Uhlmann, eds., vol. 559 of Contemporary Mathematics, Amer. Math. Soc., 
Providence, RI, 2011, pp. 1-12. 

[22] R. Dautray and J.-L. Lions, Mathematical Analysis and Numerical Methods for 
Science and Technology, Vol VI, Springer-Verlag, Berlin, 1993. 

[23] P. Elbau, O. Scherzer, and R. Schulze, Reconstruction formulas for photoacous¬ 
tic sectional imaging, Inverse Problems, 28 (2012). 045004. 

[24] D. Finch, M. Haltmeier, and Rakesh, Inversion of spherical means and the wave 
equation in even dimensions, SIAM J. Appl. Math., 68 (2007), pp. 392-412. 

[25] A. R. Fisher, A. J. Schissler, and J. C. Schotland, Photoacoustic effect for 
multiply scattered light, Phys. Rev. E, 76 (2007). 036604. 

[26] H. Gao, S. Osher, and H. Zhao, Quantitative photoacoustic tomography, in Math¬ 
ematical Modeling in Biomedical Imaging II: Optical, Ultrasound, and Opto-Acoustic 
Tomographies, H. Ammari, ed., vol. 2035 of Lecture Notes in Mathematics, Springer, 
2012, pp. 131-158. 

[27] A. Godavartya, E. M. Sevick-Muraca, and M. J. Eppstein, Three-dimensional 
fluorescence lifetime tomography, Med. Phys., 32 (2005), pp. 992-1000. 

[28] F. CoLSE, P.-L. Lions, B. Perthame, and R. Sentis, Regularity of the moments 
of the solution of a transport equation, J. Func. Anal., 76 (1988), pp. 110-125. 

[29] W. Greeberg and S. Sancaktar, Solution of the multigroup transport equations 
in IT spaces, J. Math. Phys., 17 (1976), pp. 2092-2097. 

[30] M. Haltmeier, A mollification approach for inverting the spherical mean Radon trans¬ 
form, SIAM J. Appl. Math., 71 (2011), pp. 1637-1652. 


25 



[31] M. Haltmeier, T. Schuster, and O. Scherzer, Filtered backprojection for ther¬ 
moacoustic computed tomography in spherical geometry, Math. Methods Apph Sci., 28 
(2005), pp. 1919-1937. 

[32] L. G. Henyey and J. L. Greenstein, Diffuse radiation in the galaxy, Astrophys. 
J., 90 (1941), pp. 70-83. 

[33] A. H. Hielscher, Optical tomographic imaging of small animals. Current Opinion in 
Biotechnology, 16 (2005), pp. 79-88. 

[34] Y. Hristova, Time reversal in thermoacoustic tomography - an error estimate. Inverse 
Problems, 25 (2009). 055008. 

[35] T. Kato, Perturbation Theory for Linear Operators, Springer-Verlag, Berlin, 2013. 

[36] A. Kirsch, An Introduction to the Mathematieal Theory of Inverse Problems, Springer- 
Verlag, New York, second ed., 2011. 

[37] A. Kirsch and O. Scherzer, Simultaneous reconstructions of absorption density 
and wave speed with photoacoustic measurements, SIAM J. Appl. Math., 72 (2013), 
pp. 1508-1523. 

[38] M. V. Klibanov and M. Yamamoto, Exact controllability of the time dependent 
transport equation, SIAM J. Control Optim., 46 (2007), pp. 2071-2095. 

[39] P. Kuchment, Mathematics of hybrid imaging, a brief review, in The Mathematical 
Legacy of Leon Ehrenpreis, 1. Sabadini and D. Struppa, eds., Springer-Verlag, 2012. 

[40] P. Kuchment and L. Kunyansky, Mathematics of thermoacoustic tomography, 
Euro. J. Appl. Math., 19 (2008), pp. 191-224. 

[41] -, Mathematics of thermoacoustic and photoacoustic tomography, in Handbook of 

Mathematical Methods in Imaging, O. Scherzer, ed., Springer-Verlag, 2010, pp. 817- 

866 . 

[42] A. T. N. Kumar, S. B. Raymond, A. K. Dunn, B. J. Bacskai, and D. A. 
Boas, A time domain fluorescence tomography system for small animal imaging, IEEE 
Trans. Med. Imag., 27 (2008), pp. 1152-1163. 

[43] C. Li and L. Wang, Photoacoustie tomography and sensing in biomedieine, Phys. 
Med. Biol., 54 (2009), pp. R59-R97. 

[44] A. V. Mamonov and K. Ren, Quantitative photoacoustic imaging in radiative trans¬ 
port regime, Comm. Math. Sci., 12 (2014), pp. 201-234. 

[45] M. Mokhtar-Kharroubi, Mathematical Topics in Neutron Transport Theory: New 
Aspeets, World Scientific, Singapore, 1998. 

[46] G. Y. Panasyuk, Z.-M. Wang, J. C. Schotland, and V. A. Market, Fluores¬ 
cent optical tomography with large data sets. Opt. Lett., 33 (2008), pp. 1744-1746. 


26 



[47] S. K. Patch and O. Scherzer, Photo- and thermo- acoustic imaging, Inverse Prob¬ 
lems, 23 (2007), pp. Sl-SlO. 

[48] M. S. Patterson and B. W. Pogue, Mathematical model for time-resolved and 
frequency-domain fluorescence speetroscopy in biological tissues, Appl. Opt., 33 (1994), 
pp. 1963-1974. 

[49] A. PuLKKiNEN, B. T. Cox, S. R. Arridge, J. P. Kaipio, and T. Tarvainen, A 
Bayesian approach to spectral quantitative photoacoustic tomography, Inverse Problems, 
30 (2014). 065012. 

[50] J. Qian, P. Stepanov, G. Uhlmann, and H. Zhao, An efficient Neumann-series 
based algorithm for thermoacoustic and photoacoustic tomography with variable sound 
speed, SIAM J. Imaging Sci., 4 (2011), pp. 850-883. 

[51] D. Razansky, M. Distel, C. Vinegoni, R. Ma, N. Perrimon, R. W. Koster, 
AND V. Ntziachristos, Multispectral opto-acoustic tomography of deep-seated fluo¬ 
rescent proteins in vivo. Nature Photonics, 3 (2009), pp. 412-417. 

[52] D. Razansky and V. Ntziachristos, Hybrid photoacoustic fluorescence molecular 
tomography using finite-element-based inversion, Med. Phys., 34 (2007), pp. 4293-4301. 

[53] K. Ren, Existenee and uniqueness of LA solutions to a radiative transport system, 
Preprint, (2015). 

[54] K. Ren, G. Bal, and A. H. Hielscher, Frequeney domain optieal tomography based 
on the equation of radiative transfer, SIAM J. Sci. Comput., 28 (2006), pp. 1463-1489. 

[55] -, Transport- and diffusion-based optical tomography in small domains: A compar¬ 

ative study, Applied Optics, 46 (2007), pp. 6669-6679. 

[56] K. Ren, H. Gao, and H. Zhao, A hybrid reeonstruction method for quantitative 
photoacoustic imaging, SIAM J. Imag. Sci., 6 (2013), pp. 32-55. 

[57] K. Ren and H. Zhao, Quantitative fluorescence photoacoustic tomography, SIAM J. 
Imag. Sci., 6 (2013), pp. 2024-2049. 

[58] T. Saratoon, T. Tarvainen, B. T. Cox, and S. R. Arridge, A gradient-based 
method for quantitative photoacoustic tomography using the radiative transfer equation, 
Inverse Problems, 29 (2013). 075006. 

[59] O. Scherzer, Handbook of Mathematical Methods in Imaging, Springer-Verlag, 2010. 

[60] V. Y. Soloviev, K. B. Tahir, J. McGinty, D. S. Epson, M. A. A. Neil, 
P. M. W. French, and S. R. Arridge, Fluorescence lifetime imaging by using 
time gated data acquisition. Applied Optics, 46 (2007), pp. 7384-7391. 

[61] P. Stefanov, W. Cong, and G. Wang, Modulated luminescence tomography. In¬ 
verse Problems and Imaging, 9 (2015), pp. 551-578. 


27 



[62] P. Stepanov and G. Uhlmann, Thermoacoustic tomography with variable sound 
speed, Inverse Problems, 25 (2009). 075011. 

[63] J. Tervo, On coupled Boltzmann transport equation related to radiation therapy, J. 
Math. Anal. AppL, 335 (2007), pp. 819-840. 

[64] J. Tervo and P. Kokkonen, On existence of -solutions for coupled Boltz¬ 
mann transport equation and radiation therapy treatment optimization, arXiv, (2014). 
1406.3228vl. 

[65] B. Wang, Q. Zhao, N. M. Barkey, D. L. Morse, and H. Jiang, Photoacous¬ 
tic tomography and fluorescence molecular tomography: A comparative study based on 
indocyanine green, Med. Phys., 39 (2012), pp. 2512-2517. 

[66] L. V. Wang, ed., Photoacoustic Imaging and Spectroscopy, Taylor & Francis, 2009. 

[67] -, Photoacoustic tomography, Scholarpedia, 9 (2014). 10278. 

[68] -, Photoacoustic tomography: Principles and advances, Progress in Electromagnetics 

Research, 147 (2014), pp. 1-22. 

[69] Y. Wang, K. Maslov, C. Kim, S. Hu, and L. V. Wang, Integrated photoacoustic 
and fluorescence confocal microscopy, IEEE Trans. Biomed. Eng., 57 (2010), pp. 2576- 
2578. 

[70] A. J. Welch and M. J. C. Van-Gemert, Optical-thermal Response of Laser Irra¬ 
diated Tissue, Plenum Press, New York, 1995. 

[71] B. L. Willis and C. V. M. van der Mee, Multigroup transport equations with 
nondiagonal cross-section matrices, J. Math. Phys., 27 (1986), pp. 1633-1638. 

[72] K. E. Wilson, T. Y. Wang, and J. K. Willmann, Acoustic and photoacoustic 
molecular imaging of cancer, J. Nuclear Medicine, 54 (2013), pp. 1851-1854. 

[73] R. J. Zemp, Quantitative photoacoustic tomography with multiple optieal sources. Ap¬ 
plied Optics, 49 (2010), pp. 3566-3572. 


28 



